Expectativa de vida

Laboratorio de Datos

Importacion de librerias y datos ———————————————

options(warn=-1)
library(ggplot2)
library(grid)
library(gridExtra)
library(reshape)
library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:reshape':
## 
##     rename

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

library(ggridges)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:reshape':
## 
##     rename

## The following object is masked from 'package:gridExtra':
## 
##     combine

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(gganimate)

## No renderer backend detected. gganimate will default to writing frames to separate files
## Consider installing:
## - the `gifski` package for gif output
## - the `av` package for video output
## and restarting the R session

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ lubridate 1.9.2     ✔ tibble    3.1.8
## ✔ purrr     1.0.1     ✔ tidyr     1.3.0
## ✔ readr     2.1.4

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::combine()   masks gridExtra::combine()
## ✖ tidyr::expand()    masks reshape::expand()
## ✖ dplyr::filter()    masks plotly::filter(), stats::filter()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ dplyr::rename()    masks plotly::rename(), reshape::rename()
## ✖ lubridate::stamp() masks reshape::stamp()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(class)

## 
## Attaching package: 'class'
## 
## The following object is masked from 'package:reshape':
## 
##     condense

library(tidyr)
library(caret)

## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift

library(GGally)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

datos <- read.csv("exp_de_vida.csv")

colores_continente <- c("firebrick","deepskyblue3","goldenrod1",
                        "limegreen","black")
names(colores_continente) <- c("America","Europa","Asia",
                               "Oceania","Africa")
                                   
colores_desarrollo <- c("firebrick","deepskyblue3")
names(colores_desarrollo) <- c("Developing","Developed")

Renombramos las columnas para que mantengan cierta estructura y limpiamos algunas que contienen datos imposibles o altamente sospechosos

# limpieza de columnas ---------------------------------------------------------

datos <- datos %>% 
  select(-c(BMI, percentage.expenditure, GDP, Population)) %>% 
  dplyr::rename("pais" = "Country",
         "region" = "Region",
         "ano" = "Year",
         "nivel" = "Status",
         "exp_de_vida" = "Life.expectancy",
         "mortalidad_adultos" = "Adult.Mortality",
         "muertes_infantes" = "infant.deaths",
         "alcohol" = "Alcohol",
         "hepatitis_b" = "Hepatitis.B",
         "sarampion" = "Measles",
         "muertes_menores_5" = "under.five.deaths",
         "polio" = "Polio",
         "gasto_total" = "Total.expenditure",
         "difteria" = "Diphtheria",
         "vih_sida" = "HIV.AIDS",
         "delgadez_10_a_19" = "thinness..1.19.years",
         "delgadez_5_a_9" = "thinness.5.9.years",
         "indice_icr" = "Income.composition.of.resources",
         "escolaridad" = "Schooling")

datos$muertes_infantes[datos$muertes_infantes > 800] <- NA
datos$sarampion[datos$sarampion > 900] <- NA
datos$muertes_menores_5[datos$muertes_menores_5 > 700] <- NA
datos$indice_icr[datos$indice_icr == 0] <- NA
datos$escolaridad[datos$indice_icr == 0] <- NA


datos$region <- as.character(datos$region)
datos$region[datos$region == "America del Sur"] <- "America"
datos$region[datos$region == "America del Norte"] <- "America"
datos$region[datos$region == "America Central"] <- "America"


datos$pais <- as.character(datos$pais)
datos$pais[datos$pais == "Bolivia (Plurinational State of)"] <- "Bolivia"
datos$pais[datos$pais == "Venezuela (Bolivarian Republic of)"] <- "Venezuela"
datos$pais[datos$pais == "United Kingdom of Great Britain and Northern Ireland"] <- "UK"
datos$pais[datos$pais == "United Republic of Tanzania"] <- "Tanzania"
datos$pais[datos$pais == "United States of America"] <- "USA"
datos$pais[datos$pais == "Iran (Islamic Republic of)"] <- "Iran"
datos$pais[datos$pais == "Russian Federation"] <- "Russia"
datos$pais[datos$pais == "Viet Nam"] <- "Vietnam"
datos$pais[datos$pais == "Côte d'Ivoire"] <- "Ivory Coast"
datos$pais[datos$pais == "Congo"] <- "Republic of Congo"
datos$pais[datos$pais == "Brunei Darussalam"] <- "Brunei"

# agregamos datos del GDP ------------------------------------------------------

library(tidyr)

gdp <- read.csv("gdp.csv")
gdp <- gdp[, c(1,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28)]

gdp <- gather(gdp, key = "ano", value = "Value", "X2000","X2001","X2002","X2003","X2004","X2005","X2006","X2007","X2008","X2009","X2010","X2011","X2012","X2013","X2014","X2015")
gdp$ano[gdp$ano == "X2000"] <- "2000"
gdp$ano[gdp$ano == "X2001"] <- "2001"
gdp$ano[gdp$ano == "X2002"] <- "2002"
gdp$ano[gdp$ano == "X2003"] <- "2003"
gdp$ano[gdp$ano == "X2004"] <- "2004"
gdp$ano[gdp$ano == "X2005"] <- "2005"
gdp$ano[gdp$ano == "X2006"] <- "2006"
gdp$ano[gdp$ano == "X2007"] <- "2007"
gdp$ano[gdp$ano == "X2008"] <- "2008"
gdp$ano[gdp$ano == "X2009"] <- "2009"
gdp$ano[gdp$ano == "X2010"] <- "2010"
gdp$ano[gdp$ano == "X2011"] <- "2011"
gdp$ano[gdp$ano == "X2012"] <- "2012"
gdp$ano[gdp$ano == "X2013"] <- "2013"
gdp$ano[gdp$ano == "X2014"] <- "2014"
gdp$ano[gdp$ano == "X2015"] <- "2015"

gdp$ano <- as.integer(gdp$ano)
colnames(gdp) <- c("pais","ano","gdp")

gdp$pais[gdp$pais == "United States"] <- "USA"
gdp$pais[gdp$pais == "Bahamas"] <- "Bahamas"
gdp$pais[gdp$pais == "Congo, Dem. Rep."] <- "Democratic Republic of the Congo"
gdp$pais[gdp$pais == "Congo, Rep."] <- "Republic of Congo"
gdp$pais[gdp$pais == "United Kingdom"] <- "UK"
gdp$pais[gdp$pais == "Egypt, Arab Rep."] <- "Egypt"
gdp$pais[gdp$pais == "  Gambia, The"] <- "Gambia"

datos <- left_join(datos, gdp, by=c('pais'='pais', 'ano'='ano'))

Veamos la distribucion de la variable estrella de nuestro proyecto: expectativa de vida

# distribucion de exp. de vida por region -------------------------------------- 

datos_2000 <- datos[is.element(datos$ano,"2000"),]
datos_2005 <- datos[is.element(datos$ano,"2005"),]
datos_2010 <- datos[is.element(datos$ano,"2010"),]
datos_2014 <- datos[is.element(datos$ano,"2014"),]
datos_2015 <- datos[is.element(datos$ano,"2015"),]

g2000 <- ggplot(datos_2000, aes(x = exp_de_vida, y = region, fill = region)) +
  geom_density_ridges(quantile_lines = TRUE, quantiles = 2, color="grey35") +
  theme_ridges() + 
  theme(legend.position = "none") +
  ggtitle("2000")+
  xlim(30,100)+
  ylab("")+
  xlab("exp. de vida")+
  scale_fill_manual(values = colores_continente)

g2005 <- ggplot(datos_2005, aes(x = exp_de_vida, y = region, fill = region)) +
  geom_density_ridges(quantile_lines = TRUE, quantiles = 2, color="grey35") +
  theme_ridges() + 
  theme(legend.position = "none") +
  ggtitle("2005")+
  xlim(30,100)+
  ylab("")+
  xlab("exp. de vida")+
  scale_fill_manual(values = colores_continente)

g2010 <- ggplot(datos_2010, aes(x = exp_de_vida, y = region, fill = region)) +
  geom_density_ridges(quantile_lines = TRUE, quantiles = 2, color="grey35") +
  theme_ridges() + 
  theme(legend.position = "none") +
  ggtitle("2010")+
  xlim(30,100)+
  ylab("")+
  xlab("exp. de vida")+
  scale_fill_manual(values = colores_continente)

g2015 <- ggplot(datos_2015, aes(x = exp_de_vida, y = region, fill = region)) +
  geom_density_ridges(quantile_lines = TRUE, quantiles = 2, color="grey35") +
  theme_ridges() + 
  theme(legend.position = "none") +
  ggtitle("2015")+
  xlim(30,100)+
  ylab("")+
  xlab("exp. de vida")+
  scale_fill_manual(values = colores_continente)

print(grid.arrange(g2000,g2005,g2010,g2015, nrow = 2, ncol = 2))

## Picking joint bandwidth of 2.33

## Picking joint bandwidth of 2.15

## Picking joint bandwidth of 2

## Picking joint bandwidth of 1.94

## TableGrob (2 x 2) "arrange": 4 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]
## 3 3 (2-2,1-1) arrange gtable[layout]
## 4 4 (2-2,2-2) arrange gtable[layout]

# distribucion de exp. de vida por nivel de desarrollo-------------------------- 

n2000 <- ggplot(datos_2000, aes(x = exp_de_vida, y = nivel, fill = nivel)) +
  geom_density_ridges(quantile_lines = TRUE, quantiles = 2, color="black") +
  theme_ridges() + 
  theme(legend.position = "none") +
  ggtitle("2000")+
  xlim(30,100)+
  ylab("")+
  xlab("exp. de vida")+
  scale_fill_manual(values = colores_desarrollo)

n2005 <- ggplot(datos_2005, aes(x = exp_de_vida, y = nivel, fill = nivel)) +
  geom_density_ridges(quantile_lines = TRUE, quantiles = 2, color="black") +
  theme_ridges() + 
  theme(legend.position = "none") +
  ggtitle("2005")+
  xlim(30,100)+
  ylab("")+
  xlab("exp. de vida")+
  scale_fill_manual(values = colores_desarrollo)

n2010 <- ggplot(datos_2010, aes(x = exp_de_vida, y = nivel, fill = nivel)) +
  geom_density_ridges(quantile_lines = TRUE, quantiles = 2, color="black") +
  theme_ridges() + 
  theme(legend.position = "none") +
  ggtitle("2010")+
  xlim(30,100)+
  ylab("")+
  xlab("exp. de vida")+
  scale_fill_manual(values = colores_desarrollo)

n2015 <- ggplot(datos_2015, aes(x = exp_de_vida, y = nivel, fill = nivel)) +
  geom_density_ridges(quantile_lines = TRUE, quantiles = 2, color="black") +
  theme_ridges() + 
  theme(legend.position = "none") +
  ggtitle("2015")+
  xlim(30,100)+
  ylab("")+
  xlab("exp. de vida")+
  scale_fill_manual(values = colores_desarrollo)

grid.arrange(n2000,n2005,n2010,n2015, nrow = 2, ncol = 2)

## Picking joint bandwidth of 2.14

## Picking joint bandwidth of 2.18

## Picking joint bandwidth of 2.23

## Picking joint bandwidth of 1.93

# distribucion de nivel de desarrollo por continente --------------------------- 

ggplot(datos_2015, aes(x = region, fill = nivel)) +
  geom_bar(position="fill")+
  ggtitle("Proporción del nivel de desarrollo por continente
")+
  ylab("proporcion")+
  xlab("region")+
  scale_fill_manual(values = colores_desarrollo)+
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))

# el tamano de las burbujas depende del pbi (mas grande, mas pbi) --------------
# el tamano de las burbujas depende del pbi (mas grande, mas pbi) --------------
# el tamano de las burbujas depende del pbi (mas grande, mas pbi) --------------
# el tamano de las burbujas depende del pbi (mas grande, mas pbi) --------------
# el tamano de las burbujas depende del pbi (mas grande, mas pbi) --------------
# el tamano de las burbujas depende del pbi (mas grande, mas pbi) --------------
# el tamano de las burbujas depende del pbi (mas grande, mas pbi) --------------
# el tamano de las burbujas depende del pbi (mas grande, mas pbi) --------------


r3 <- plot_ly(datos, x = ~`indice_icr`, y = ~`exp_de_vida`, ylab="exp. de vida", xlab="gasto total",
              frame=~ano,size= ~`gdp`,color=~`region`,colors=colores_continente,mode="markers",marker = list(symbol = 'circle', sizemode = 'diameter', line = list(width = 2, color = '#FFFFFF'), opacity=0.75)) %>% animation_opts(1000, redraw = FALSE)

r3

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

# Gasto total vs exp. de vida  -------------------------------------------------

v1 <- ggplot(datos_2000,aes(gasto_total,exp_de_vida,color = region))+
  geom_point()+
  ggtitle("2000")+
  xlab("Gasto total")+
  ylab("Expectativa de vida")+
  xlim(0,18)+
  ylim(35,90)+
  scale_color_manual(values= colores_continente)
  
v2 <- ggplot(datos_2005,aes(gasto_total,exp_de_vida,color = region))+
  geom_point()+
  ggtitle("2005")+
  xlab("Gasto total")+
  ylab("Expectativa de vida")+
  xlim(0,18)+
  ylim(35,90)+
  scale_color_manual(values= colores_continente)

v3 <- ggplot(datos_2010,aes(gasto_total,exp_de_vida,color = region))+
  geom_point()+
  ggtitle("2010")+
  xlab("Gasto total")+
  ylab("Expectativa de vida")+
  xlim(0,18)+
  ylim(35,90)+
  scale_color_manual(values= colores_continente)

v4 <- ggplot(datos_2014,aes(gasto_total,exp_de_vida,color = region))+
  geom_point()+
  ggtitle("2014")+
  xlab("Gasto total")+
  ylab("Expectativa de vida")+
  xlim(0,18)+
  ylim(35,90)+
  scale_color_manual(values= colores_continente)


grid.arrange(v1,v2,v3,v4,nrow=2,ncol=2)

Tambien veamoslo en los mapas

mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2000, by="pais")
mapdata2000<-mapdata %>% filter(!is.na(mapdata$exp_de_vida))

mapa2000<-ggplot(mapdata2000, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2000")+
  scale_fill_gradient(name = "exp. de vida", low = "gold1", high =  "navyblue", limits = c(40,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())



mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- full_join(mapdata, datos_2010, by="pais")
mapdata2005<-mapdata %>% filter(!is.na(mapdata$exp_de_vida))

mapa2005<-ggplot(mapdata2005, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2005")+
  scale_fill_gradient(name = "exp. de vida", low = "gold1", high =  "navyblue", limits = c(40,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- full_join(mapdata, datos_2010, by="pais")
mapdata2010<-mapdata %>% filter(!is.na(mapdata$exp_de_vida))

mapa2010<-ggplot(mapdata2010, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2010")+
  scale_fill_gradient(name = "exp. de vida", low = "gold1", high =  "navyblue", limits = c(40,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- full_join(mapdata, datos_2010, by="pais")
mapdata2015<-mapdata %>% filter(!is.na(mapdata$exp_de_vida))

mapa2015<-ggplot(mapdata2015, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2015")+
  scale_fill_gradient(name = "exp. de vida", low = "gold1", high =  "navyblue", limits = c(40,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())




grid.arrange(mapa2000,mapa2005,mapa2010,mapa2015, ncol=2,nrow=2)

# EUROPA --------

# graficos ------

mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2000, by="pais")
mapdata2000<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Europa"))

mapa2000<-ggplot(mapdata2000, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2000")+
  scale_fill_gradient(name = "exp. de vida", low = "cadetblue1", high =  "darkblue", limits = c(65,90), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2005, by="pais")
mapdata2005<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Europa"))

mapa2005<-ggplot(mapdata2005, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2005")+
  scale_fill_gradient(name = "exp. de vida", low = "cadetblue1", high =  "darkblue", limits = c(65,90), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2010, by="pais")
mapdata2010<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Europa"))

mapa2010<-ggplot(mapdata2010, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2010")+
  scale_fill_gradient(name = "exp. de vida", low = "cadetblue1", high =  "darkblue", limits = c(65,90), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())



mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2015, by="pais")
mapdata2015<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Europa") & !is.element(mapdata$pais, "Bahamas"))

mapa2015<-ggplot(mapdata2015, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2015")+
  scale_fill_gradient(name = "exp. de vida", low = "cadetblue1", high =  "darkblue", limits = c(65,90), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())

grid.arrange(mapa2000,mapa2005,mapa2010,mapa2015, nrow=2, ncol=2)

# ASIA ----------

# graficos ------

mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2000, by="pais")
mapdata2000<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Asia"))

mapa2000<-ggplot(mapdata2000, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2000")+
  scale_fill_gradient(name = "exp. de vida", low = "khaki1", high =  "saddlebrown", limits = c(60,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2005, by="pais")
mapdata2005<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Asia"))

mapa2005<-ggplot(mapdata2005, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2005")+
  scale_fill_gradient(name = "exp. de vida", low = "khaki1", high =  "saddlebrown", limits = c(60,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2010, by="pais")
mapdata2010<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Asia"))

mapa2010<-ggplot(mapdata2010, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2010")+
  scale_fill_gradient(name = "exp. de vida", low = "khaki1", high =  "saddlebrown", limits = c(60,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())



mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2015, by="pais")
mapdata2015<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Asia"))

mapa2015<-ggplot(mapdata2015, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2015")+
  scale_fill_gradient(name = "exp. de vida", low = "khaki1", high =  "saddlebrown", limits = c(60,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())

grid.arrange(mapa2000,mapa2005,mapa2010,mapa2015, nrow=2, ncol=2)

# AFRICA ---------


# graficos ------

mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2000, by="pais")
mapdata2000<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Africa"))

mapa2000<-ggplot(mapdata2000, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "white")+
  ggtitle("Exp. de vida en el 2000")+
  scale_fill_gradient(name = "exp. de vida", low = "gray80", high =  "black", limits = c(40,80), na.value = "lightyellow")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2005, by="pais")
mapdata2005<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Africa"))

mapa2005<-ggplot(mapdata2005, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "white")+
  ggtitle("Exp. de vida en el 2005")+
  scale_fill_gradient(name = "exp. de vida", low = "gray80", high =  "black", limits = c(40,80), na.value = "lightyellow")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2010, by="pais")
mapdata2010<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Africa"))

mapa2010<-ggplot(mapdata2010, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "white")+
  ggtitle("Exp. de vida en el 2010")+
  scale_fill_gradient(name = "exp. de vida", low = "gray80", high =  "black", limits = c(40,80), na.value = "lightyellow")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())



mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2015, by="pais")
mapdata2015<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "Africa"))

mapa2015<-ggplot(mapdata2015, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "white")+
  ggtitle("Exp. de vida en el 2015")+
  scale_fill_gradient(name = "exp. de vida", low = "gray80", high =  "black", limits = c(40,80), na.value = "lightyellow")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())

grid.arrange(mapa2000,mapa2005,mapa2010,mapa2015, nrow=2, ncol=2)

# AMERICA -------
  


# graficos ------

mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2000, by="pais")
mapdata2000<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "America"))

mapa2000<-ggplot(mapdata2000, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2000")+
  scale_fill_gradient(name = "exp. de vida", low = "white", high =  "darkred", limits = c(60,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2005, by="pais")
mapdata2005<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "America"))

mapa2005<-ggplot(mapdata2005, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2005")+
  scale_fill_gradient(name = "exp. de vida", low = "white", high =  "darkred", limits = c(60,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())


mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2010, by="pais")
mapdata2010<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "America"))

mapa2010<-ggplot(mapdata2010, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2010")+
  scale_fill_gradient(name = "exp. de vida", low = "white", high =  "darkred", limits = c(60,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())



mapdata <- map_data("world") ##ggplot2
colnames(mapdata)[5] <- "pais"
mapdata <- left_join(mapdata, datos_2015, by="pais")
mapdata2015<-mapdata %>% filter(!is.na(mapdata$exp_de_vida) & is.element(mapdata$region, "America"))

mapa2015<-ggplot(mapdata2015, aes( x = long, y = lat, group=group)) +
  geom_polygon(aes(fill = exp_de_vida), color = "grey10")+
  ggtitle("Exp. de vida en el 2015")+
  scale_fill_gradient(name = "exp. de vida", low = "white", high =  "darkred", limits = c(60,85), na.value = "grey")+
  theme(axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.y=element_blank(),
        axis.title.x=element_blank(),
        rect = element_blank())

grid.arrange(mapa2000,mapa2005,mapa2010,mapa2015, nrow=2, ncol=2)

datos_exp <- datos[!is.na(datos$exp_de_vida),]
promedios_exp <- datos_exp%>% group_by(region,ano)%>%summarise(mean_val=mean(exp_de_vida))

## `summarise()` has grouped output by 'region'. You can override using the
## `.groups` argument.

ggplot(promedios_exp, aes(ano,mean_val,color=region))+
  geom_line(size=1)+
  xlab("Ano")+
  ylab("Promedio de exp. de vida")+
  scale_color_manual(values = colores_continente)

america <- datos[is.element(datos$region, "America"),]
paises_america <- unique(america$pais)

for(i in paises_america){
  if(i=="Haiti"){
    color_if_haiti <- "firebrick"
  }else{
    color_if_haiti <- "black"
  }
  
  
  data <- america[america$pais== i,]
  print(ggplot(data, aes(ano,exp_de_vida,color=pais))+
    geom_line(size=1)+
    xlab("Ano")+
    ylab("Promedio de exp. de vida")+
    ggtitle(i)+
    geom_vline(xintercept = 2010)+
    scale_color_manual(values=color_if_haiti)
  )
  
}

# Modelo America ------------------------------------------------------------------------------

america <- datos[is.element(datos$region, "America"),]
america$gasto_total[is.na(america$gasto_total)] <- mean(america$gasto_total[-1])
america$alcohol[is.na(america$alcohol)] <- mean(america$alcohol[-c(1)])
america$gdp[is.na(america$gdp)] <- mean(america$gdp[-c(1)])

# correlacion entre variables (america)

ggpairs(america[,c(5,18,19,20)], title="correlogram with ggpairs()") 

set.seed(73) # Fijamos semilla para que siempre retorne el mismo resultado


ggplot(america,aes(escolaridad,exp_de_vida))+
  geom_point()+
  xlab("Escolaridad")+
  ylab("Expectativa de vida")

ggplot(america,aes(indice_icr,exp_de_vida))+
  geom_point()+
  xlab("Indice ICR")+
  ylab("Expectativa de vida")

ggplot(america,aes(gdp,exp_de_vida))+
  geom_point()+
  xlab("GDP")+
  ylab("Expectativa de vida")

# Modelo exp. vida y gdp -------------------------------------------------------

america_modelo <- america[!is.na(america$exp_de_vida) & !is.na(america$gdp),]

Y <- america_modelo$exp_de_vida
L1 <- america_modelo$gdp

Eval <- function(mu, alfa) {
  salida<-mean((Y-mu-alfa*L1)^2)
  return(t(salida))
}
Eval(10,2)

##            [,1]
## [1,] 1039979996

facred<-0.9999 # factor de reduccion de la ventana
facred.acu<-1 # factor de reduccion acumulado
toler<-1/1e5 # umbral de tolerancia
# rangos
rango.mu<- 50 # rango inicial de mu
rango.alfa<- 25 # rango inicial de alfa
# parametros iniciales
mu<- 20 # valor inicial de mu
alfa<- 10 # valor inicial de alfa
# parametros mejores
mejor.mu<-mu
mejor.alfa<-alfa
mejor.eval<-Eval(mu,alfa) # el mejor valor
mejores<-matrix(c(mejor.eval,mu,alfa,1),1,5)
k<-0 # indice de iteracion
actu<-0 # indice de actualizacion


while (facred.acu>toler)
{
  k<-k+1
  # Genero nuevos valores aleatorios
  mu<-runif(1,mejor.mu-rango.mu*facred.acu,mejor.mu+rango.mu*facred.acu)
  alfa<-runif(1,mejor.alfa-rango.alfa*facred.acu,mejor.alfa+rango.alfa*facred.acu)
  # Evaluacion de los nuevos valores
  valor<-Eval(mu,alfa)
  if (valor<mejor.eval) # SI encuentro algo mejor -> Actualizacion
  {
    actu<-actu+1
    mejor.eval<-valor
    mejor.mu<-mu
    mejor.alfa<-alfa
    mejores<-rbind(mejores,c(mejor.eval,mu,alfa,k))
  }
  else # SI NO encuentro algo mejor -> Reduzco rango de busqueda
  {
    facred.acu<-facred.acu*facred
  }
}

c(mu,alfa)

## [1] 25.184196111  0.002256684

ajusA1_1<-lm(Y~L1,data=america_modelo) # modelo lineal multiple

predichos<-predict(ajusA1_1)
mean(abs(Y-predichos))

## [1] 2.883065

pmaeA1_1<-mean(abs(Y-predichos))/mean(Y)
pmaeA1_1

## [1] 0.03947337

n<-length(Y)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Y~L1,data=america_modelo[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=america_modelo[i,])
}
mean(abs(Y-predichos.oos))

## [1] 3.926512

pmaeA1_1.cv<-mean(abs(Y-predichos.oos))/mean(Y)
pmaeA1_1.cv

## [1] 0.05375969

# Modelo exp. vida y escolaridad------------------------------------------------

america_modelo <- america[!is.na(america$exp_de_vida) & !is.na(america$escolaridad),]

Y <- america_modelo$exp_de_vida
L1 <- america_modelo$escolaridad

Eval <- function(mu, alfa) {
  salida<-mean((Y-mu-alfa*L1)^2)
  return(t(salida))
}
Eval(10,2)

##          [,1]
## [1,] 1439.139

facred<-0.9999 # factor de reduccion de la ventana
facred.acu<-1 # factor de reduccion acumulado
toler<-1/1e5 # umbral de tolerancia
# rangos
rango.mu<- 50 # rango inicial de mu
rango.alfa<- 25 # rango inicial de alfa
# parametros iniciales
mu<- 20 # valor inicial de mu
alfa<- 10 # valor inicial de alfa
# parametros mejores
mejor.mu<-mu
mejor.alfa<-alfa
mejor.eval<-Eval(mu,alfa) # el mejor valor
mejores<-matrix(c(mejor.eval,mu,alfa,1),1,5)
k<-0 # indice de iteracion
actu<-0 # indice de actualizacion


while (facred.acu>toler)
{
  k<-k+1
  # Genero nuevos valores aleatorios
  mu<-runif(1,mejor.mu-rango.mu*facred.acu,mejor.mu+rango.mu*facred.acu)
  alfa<-runif(1,mejor.alfa-rango.alfa*facred.acu,mejor.alfa+rango.alfa*facred.acu)
  # Evaluacion de los nuevos valores
  valor<-Eval(mu,alfa)
  if (valor<mejor.eval) # SI encuentro algo mejor -> Actualizacion
  {
    actu<-actu+1
    mejor.eval<-valor
    mejor.mu<-mu
    mejor.alfa<-alfa
    mejores<-rbind(mejores,c(mejor.eval,mu,alfa,k))
  }
  else # SI NO encuentro algo mejor -> Reduzco rango de busqueda
  {
    facred.acu<-facred.acu*facred
  }
}

c(mu,alfa)

## [1] 60.6258465  0.9801511

ajusA1_2<-lm(Y~L1,data=america_modelo) # modelo lineal multiple

predichos<-predict(ajusA1_2)
mean(abs(Y-predichos))

## [1] 2.826704

pmaeA1_2<-mean(abs(Y-predichos))/mean(Y)
pmaeA1_2

## [1] 0.03867492

n<-length(Y)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Y~L1,data=america_modelo[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=america_modelo[i,])
}
mean(abs(Y-predichos.oos))

## [1] 3.214812

pmaeA1_2.cv<-mean(abs(Y-predichos.oos))/mean(Y)
pmaeA1_2.cv

## [1] 0.04398501

# Modelo exp. vida e icr--------------------------------------------------------

america_modelo <- america[!is.na(america$exp_de_vida) & !is.na(america$indice_icr),]

Y <- america_modelo$exp_de_vida
L1 <- america_modelo$indice_icr

Eval <- function(mu, alfa) {
  salida<-mean((Y-mu-alfa*L1)^2)
  return(t(salida))
}
Eval(10,2)

##          [,1]
## [1,] 3829.775

facred<-0.9999 # factor de reduccion de la ventana
facred.acu<-1 # factor de reduccion acumulado
toler<-1/1e5 # umbral de tolerancia
# rangos
rango.mu<- 50 # rango inicial de mu
rango.alfa<- 25 # rango inicial de alfa
# parametros iniciales
mu<- 20 # valor inicial de mu
alfa<- 10 # valor inicial de alfa
# parametros mejores
mejor.mu<-mu
mejor.alfa<-alfa
mejor.eval<-Eval(mu,alfa) # el mejor valor
mejores<-matrix(c(mejor.eval,mu,alfa,1),1,5)
k<-0 # indice de iteracion
actu<-0 # indice de actualizacion


while (facred.acu>toler)
{
  k<-k+1
  # Genero nuevos valores aleatorios
  mu<-runif(1,mejor.mu-rango.mu*facred.acu,mejor.mu+rango.mu*facred.acu)
  alfa<-runif(1,mejor.alfa-rango.alfa*facred.acu,mejor.alfa+rango.alfa*facred.acu)
  # Evaluacion de los nuevos valores
  valor<-Eval(mu,alfa)
  if (valor<mejor.eval) # SI encuentro algo mejor -> Actualizacion
  {
    actu<-actu+1
    mejor.eval<-valor
    mejor.mu<-mu
    mejor.alfa<-alfa
    mejores<-rbind(mejores,c(mejor.eval,mu,alfa,k))
  }
  else # SI NO encuentro algo mejor -> Reduzco rango de busqueda
  {
    facred.acu<-facred.acu*facred
  }
}

c(mu,alfa)

## [1] 44.21452 41.43493

ajusA1_3<-lm(Y~L1,data=america_modelo) # modelo lineal multiple

predichos<-predict(ajusA1_3)
mean(abs(Y-predichos))

## [1] 2.087726

pmaeA1_3<-mean(abs(Y-predichos))/mean(Y)
pmaeA1_3

## [1] 0.02855348

n<-length(Y)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Y~L1,data=america_modelo[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=america_modelo[i,])
}
mean(abs(Y-predichos.oos))

## [1] 4.341871

pmaeA1_3.cv<-mean(abs(Y-predichos.oos))/mean(Y)
pmaeA1_3.cv

## [1] 0.05938305

# Modelo exp. vida, icr y escolaridad ------------------------------------------

america_modelo <- america[!is.na(america$exp_de_vida) & !is.na(america$indice_icr) & !is.na(america$escolaridad),]

Y <- america_modelo$exp_de_vida
L1 <- america_modelo$indice_icr
L2 <- america_modelo$escolaridad

Eval <- function(mu, alfa, beta) {
  salida<-mean((Y-mu-alfa*L1-beta*L2)^2)
  return(t(salida))
}
Eval(10,2,1)

##          [,1]
## [1,] 2402.706

facred<-0.9999 # factor de reduccion de la ventana
facred.acu<-1 # factor de reduccion acumulado
toler<-1/1e5 # umbral de tolerancia
# rangos
rango.mu<- 50 # rango inicial de mu
rango.alfa<- 25 # rango inicial de alfa
rango.beta<- 25 # rango inicial de beta
# parametros iniciales
mu<- 20 # valor inicial de mu
alfa<- 10 # valor inicial de alfa
beta<- 10 # valor inicial de beta
# parametros mejores
mejor.mu<-mu
mejor.alfa<-alfa
mejor.beta<-beta
mejor.eval<-Eval(mu,alfa,beta) # el mejor valor
mejores<-matrix(c(mejor.eval,mu,alfa,beta,1),1,5)
k<-0 # indice de iteracion
actu<-0 # indice de actualizacion


while (facred.acu>toler)
{
  k<-k+1
  # Genero nuevos valores aleatorios
  mu<-runif(1,mejor.mu-rango.mu*facred.acu,mejor.mu+rango.mu*facred.acu)
  alfa<-runif(1,mejor.alfa-rango.alfa*facred.acu,mejor.alfa+rango.alfa*facred.acu)
  beta<-runif(1,mejor.beta-rango.beta*facred.acu,mejor.beta+rango.beta*facred.acu)
  # Evaluacion de los nuevos valores
  valor<-Eval(mu,alfa,beta)
  if (valor<mejor.eval) # SI encuentro algo mejor -> Actualizacion
  {
    actu<-actu+1
    mejor.eval<-valor
    mejor.mu<-mu
    mejor.alfa<-alfa
    mejor.beta<-beta
    mejores<-rbind(mejores,c(mejor.eval,mu,alfa,beta,k))
  }
  else # SI NO encuentro algo mejor -> Reduzco rango de busqueda
  {
    facred.acu<-facred.acu*facred
  }
}

c(mu,alfa,beta)

## [1] 44.2573072 43.6603281 -0.1239508

ajusA2_1<-lm(Y~L1+L2,data=america_modelo) # modelo lineal multiple

predichos<-predict(ajusA2_1)
mean(abs(Y-predichos))

## [1] 2.069614

pmaeA2_1<-mean(abs(Y-predichos))/mean(Y)
pmaeA2_1

## [1] 0.02830577

n<-length(Y)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Y~L1+L2,data=america_modelo[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=america_modelo[i,])
}
mean(abs(Y-predichos.oos))

## [1] 4.383204

pmaeA2_1.cv<-mean(abs(Y-predichos.oos))/mean(Y)
pmaeA2_1.cv

## [1] 0.05994835

# Modelo exp. vida, icr y gdp --------------------------------------------------

america_modelo <- america[!is.na(america$exp_de_vida) & !is.na(america$indice_icr) & !is.na(america$gdp),]

Y <- america_modelo$exp_de_vida
L1 <- america_modelo$indice_icr
L2 <- america_modelo$gdp

Eval <- function(mu, alfa, beta) {
  salida<-mean((Y-mu-alfa*L1-beta*L2)^2)
  return(t(salida))
}
Eval(10,2,1)

##           [,1]
## [1,] 192849030

facred<-0.9999 # factor de reduccion de la ventana
facred.acu<-1 # factor de reduccion acumulado
toler<-1/1e5 # umbral de tolerancia
# rangos
rango.mu<- 50 # rango inicial de mu
rango.alfa<- 25 # rango inicial de alfa
rango.beta<- 25 # rango inicial de beta
# parametros iniciales
mu<- 20 # valor inicial de mu
alfa<- 10 # valor inicial de alfa
beta<- 10 # valor inicial de beta
# parametros mejores
mejor.mu<-mu
mejor.alfa<-alfa
mejor.beta<-beta
mejor.eval<-Eval(mu,alfa,beta) # el mejor valor
mejores<-matrix(c(mejor.eval,mu,alfa,beta,1),1,5)
k<-0 # indice de iteracion
actu<-0 # indice de actualizacion


while (facred.acu>toler)
{
  k<-k+1
  # Genero nuevos valores aleatorios
  mu<-runif(1,mejor.mu-rango.mu*facred.acu,mejor.mu+rango.mu*facred.acu)
  alfa<-runif(1,mejor.alfa-rango.alfa*facred.acu,mejor.alfa+rango.alfa*facred.acu)
  beta<-runif(1,mejor.beta-rango.beta*facred.acu,mejor.beta+rango.beta*facred.acu)
  # Evaluacion de los nuevos valores
  valor<-Eval(mu,alfa,beta)
  if (valor<mejor.eval) # SI encuentro algo mejor -> Actualizacion
  {
    actu<-actu+1
    mejor.eval<-valor
    mejor.mu<-mu
    mejor.alfa<-alfa
    mejor.beta<-beta
    mejores<-rbind(mejores,c(mejor.eval,mu,alfa,beta,k))
  }
  else # SI NO encuentro algo mejor -> Reduzco rango de busqueda
  {
    facred.acu<-facred.acu*facred
  }
}

c(mu,alfa,beta)

## [1] 33.982983589 32.262766476  0.001105168

ajusA2_2<-lm(Y~L1+L2,data=america_modelo) # modelo lineal multiple

predichos<-predict(ajusA2_2)
mean(abs(Y-predichos))

## [1] 2.07257

pmaeA2_2<-mean(abs(Y-predichos))/mean(Y)
pmaeA2_2

## [1] 0.02843491

n<-length(Y)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Y~L1+L2,data=america_modelo[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=america_modelo[i,])
}
mean(abs(Y-predichos.oos))

## [1] 4.320345

pmaeA2_2.cv<-mean(abs(Y-predichos.oos))/mean(Y)
pmaeA2_2.cv

## [1] 0.05927356

# Modelo exp. vida, escolaridad y gdp --------------------------------------------------

america_modelo <- america[!is.na(america$exp_de_vida) & !is.na(america$escolaridad) & !is.na(america$gdp),]

Y <- america_modelo$exp_de_vida
L1 <- america_modelo$escolaridad
L2 <- america_modelo$gdp

Eval <- function(mu, alfa, beta) {
  salida<-mean((Y-mu-alfa*L1-beta*L2)^2)
  return(t(salida))
}
Eval(10,2,1)

##           [,1]
## [1,] 191051331

facred<-0.9999 # factor de reduccion de la ventana
facred.acu<-1 # factor de reduccion acumulado
toler<-1/1e5 # umbral de tolerancia
# rangos
rango.mu<- 50 # rango inicial de mu
rango.alfa<- 25 # rango inicial de alfa
rango.beta<- 25 # rango inicial de beta
# parametros iniciales
mu<- 20 # valor inicial de mu
alfa<- 10 # valor inicial de alfa
beta<- 10 # valor inicial de beta
# parametros mejores
mejor.mu<-mu
mejor.alfa<-alfa
mejor.beta<-beta
mejor.eval<-Eval(mu,alfa,beta) # el mejor valor
mejores<-matrix(c(mejor.eval,mu,alfa,beta,1),1,5)
k<-0 # indice de iteracion
actu<-0 # indice de actualizacion


while (facred.acu>toler)
{
  k<-k+1
  # Genero nuevos valores aleatorios
  mu<-runif(1,mejor.mu-rango.mu*facred.acu,mejor.mu+rango.mu*facred.acu)
  alfa<-runif(1,mejor.alfa-rango.alfa*facred.acu,mejor.alfa+rango.alfa*facred.acu)
  beta<-runif(1,mejor.beta-rango.beta*facred.acu,mejor.beta+rango.beta*facred.acu)
  # Evaluacion de los nuevos valores
  valor<-Eval(mu,alfa,beta)
  if (valor<mejor.eval) # SI encuentro algo mejor -> Actualizacion
  {
    actu<-actu+1
    mejor.eval<-valor
    mejor.mu<-mu
    mejor.alfa<-alfa
    mejor.beta<-beta
    mejores<-rbind(mejores,c(mejor.eval,mu,alfa,beta,k))
  }
  else # SI NO encuentro algo mejor -> Reduzco rango de busqueda
  {
    facred.acu<-facred.acu*facred
  }
}

c(mu,alfa,beta)

## [1] -50.535691528   9.684760843  -0.000301152

ajusA2_3<-lm(Y~L1+L2,data=america_modelo) # modelo lineal multiple

predichos<-predict(ajusA2_3)
mean(abs(Y-predichos))

## [1] 2.744638

pmaeA2_3<-mean(abs(Y-predichos))/mean(Y)
pmaeA2_3

## [1] 0.03767086

n<-length(Y)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Y~L1+L2,data=america_modelo[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=america_modelo[i,])
}
mean(abs(Y-predichos.oos))

## [1] 4.611299

pmaeA2_3.cv<-mean(abs(Y-predichos.oos))/mean(Y)
pmaeA2_3.cv

## [1] 0.06329126

# Modelo exp. vida, escolaridad, icr y gdp -------------------------------------

america_modelo <- america[!is.na(america$exp_de_vida) & !is.na(america$indice_icr) &
                            !is.na(america$escolaridad) & !is.na(america$gdp),]

Y <- america_modelo$exp_de_vida
L1 <- america_modelo$indice_icr
L2 <- america_modelo$escolaridad
L3 <- america_modelo$gdp

mu<- 35 # valor inicial de mu
a<- 11 # valor inicial de a
b<- 10 # valor inicial de b
c<- 8

Eval_arg <- function(mu,a,b,c) {
  salida<-mean((Y-mu-a*L1-b*L2-c*L3)^2)
  return(t(salida))
}

Eval_arg(mu,a,b,c)

##             [,1]
## [1,] 12452691555

facred<-0.9999 # factor de reduccion de la ventana
facred.acu<-1 # factor de reduccion acumulado
toler<-1/1e3 # umbral de tolerancia
# rangos
rango.mu<- 50 # rango inicial de mu
rango.a<- 25 # rango inicial de alfa
rango.b<- 25 # rango inicial de beta
rango.c<- 25

facred.acu<-1

mejor.mu<-mu
mejor.a<-a
mejor.b<-b
mejor.c<-c
mejor.eval<-Eval_arg(mu,a,b,c) # el mejor valor
mejores<-matrix(c(mejor.eval,mu,a,b,c,1),1,5)

k<-0 # indice de iteracion
actu<-0 # indice de actualizacion

while (facred.acu>toler)
{
  k<-k+1
  # Genero nuevos valores aleatorios
  mu<-runif(1,mejor.mu-rango.mu*facred.acu,mejor.mu+rango.mu*facred.acu)
  a<-runif(1,mejor.a-rango.a*facred.acu,mejor.a+rango.a*facred.acu)
  b<-runif(1,mejor.b-rango.b*facred.acu,mejor.b+rango.b*facred.acu)
  c<-runif(1,mejor.c-rango.c*facred.acu,mejor.c+rango.c*facred.acu)
  
  # Evaluacion de los nuevos valores
  valor<-Eval_arg(mu,a,b,c)
  if (valor<mejor.eval) # SI encuentro algo mejor -> Actualizacion
  {
    actu<-actu+1
    mejor.eval<-valor
    mejor.mu<-mu
    mejor.a<-a
    mejor.b<-b
    mejor.c<-c
    mejores<-rbind(mejores,c(mejor.eval,mu,a,b,c,k))
  }
  else # SI NO encuentro algo mejor -> Reduzco rango de busqueda
  {
    facred.acu<-facred.acu*facred
  }
}

ajusA3<-lm(Y~L1+L2+L3,data=america_modelo) # modelo lineal multiple
coe<-coef(ajusA3) # coeficientes
coe

##   (Intercept)            L1            L2            L3 
## 37.1720895551 66.8000479992 -0.6775598509 -0.0001624293

c(mu,a,b,c)

## [1] 47.06674074 13.50611557  3.97565677 -0.02558169

predichos<-predict(ajusA3)
pmaeA3 <- mean(abs(Y-predichos))/mean(Y)
pmaeA3

## [1] 0.0263972

n<-length(Y)
predichos.oos<-rep(NA,n) # predichos out of sample
for (i in 1:n)
{
  ajus.cv<-lm(Y~L1+L2+L3,data=america_modelo[-i,])
  predichos.oos[i]<-predict(ajus.cv,newdata=america_modelo[i,])
}
mean(abs(Y-predichos.oos))

## [1] 4.440534

pmaeA3.cv<-mean(abs(Y-predichos.oos))/mean(Y)
pmaeA3.cv

## [1] 0.0609225

# PMAes ingenuos y PMAES CV ----------------------------------------------------


PMAESA<-matrix(c(pmaeA1_1,pmaeA1_1.cv,"GDP",
                 pmaeA1_2,pmaeA1_2.cv,"esc.",
                 pmaeA1_3,pmaeA1_3.cv,"ICR",
                 pmaeA2_1,pmaeA2_1.cv,"esc. e ICR",
                 pmaeA2_2,pmaeA2_2.cv,"GDP e ICR",
                 pmaeA2_3,pmaeA2_3.cv,"GDP y esc.",
                 pmaeA3,pmaeA3.cv,"GDP, esc., ICR"
                 ),7,3,byrow = T)


PMAESA <- as.data.frame(PMAESA)
PMAESA$V2 <- as.numeric(PMAESA$V2)
PMAESA<- PMAESA %>% arrange(desc(V2))

barplot(height = PMAESA$V2, names = PMAESA$V3, las=2, ylim=c(0,0.07), col="firebrick4", border = "firebrick4", cex.names=0.7, main="PMAEs segun el modelo utilizado (ingenuo)")

PMAESA <- as.data.frame(PMAESA)
PMAESA$V1 <- as.numeric(PMAESA$V1)
PMAESA<- PMAESA %>% arrange(desc(V1))

barplot(height = PMAESA$V1, names = PMAESA$V3, las=2, ylim=c(0,0.07), col="firebrick", border = "firebrick", cex.names=0.7, main="PMAEs segun el modelo utilizado (CV)")

# knn de exp. de vida vs gdp, predice nivel de desarrollo ----------------------

# knn de exp. de vida vs gdp, predice nivel de desarrollo ----------------------

ggplot(datos_2010, aes(nivel, exp_de_vida, fill=nivel))+
  geom_boxplot()+
  ylab("exp. de vida")+
  scale_fill_manual(values = colores_desarrollo)

ggplot(datos_2010, aes(nivel, gdp, fill=nivel))+
  geom_boxplot()+
  scale_fill_manual(values = colores_desarrollo)

y2000<-ggplot(datos_2000, aes(exp_de_vida, gdp, color=nivel))+
  geom_point()+
  ggtitle("Exp. de vida vs GDP, año 2000")+
  xlab("exp. de vida")+
  xlim(35,90)+
  scale_color_manual(values = colores_desarrollo)

y2005<-ggplot(datos_2005, aes(exp_de_vida, gdp, color=nivel))+
  geom_point()+
  ggtitle("Exp. de vida vs GDP, año 2005")+
  xlab("exp. de vida")+
  xlim(35,90)+
  scale_color_manual(values = colores_desarrollo)

y2010<-ggplot(datos_2010, aes(exp_de_vida, gdp, color=nivel))+
  geom_point()+
  ggtitle("Exp. de vida vs GDP, año 2010")+
  xlab("exp. de vida")+
  xlim(35,90)+
  scale_color_manual(values = colores_desarrollo)

y2015<-ggplot(datos_2015, aes(exp_de_vida, gdp, color=nivel))+
  geom_point()+
  ggtitle("Exp. de vida vs GDP, año 2015")+
  xlab("exp. de vida")+
  xlim(35,90)+
  scale_color_manual(values = colores_desarrollo)

grid.arrange(y2000,y2005,y2010,y2015, ncol=2, nrow=2)

datos_knn <- datos[!is.na(datos$exp_de_vida) & !is.na(datos$gdp), c(3,4,5,20)]

datos_knn["exp_de_vida"] <- (datos_knn$exp_de_vida - mean(datos_knn$exp_de_vida))/ sd(datos_knn$exp_de_vida)
datos_knn["gdp"] <- (datos_knn$gdp - mean(datos_knn$gdp))/ sd(datos_knn$gdp)


datos_knn_2000 <- datos_knn[is.element(datos_knn$ano,"2000"),]
datos_knn_2005 <- datos_knn[is.element(datos_knn$ano,"2005"),]
datos_knn_2010 <- datos_knn[is.element(datos_knn$ano,"2010"),]
datos_knn_2015 <- datos_knn[is.element(datos_knn$ano,"2015"),]


z2000<-ggplot(datos_knn_2000, aes(exp_de_vida, gdp, color=nivel))+
  geom_point()+
  ggtitle("Exp. de vida vs GDP, año 2000")+
  xlim(-3.2,2.2)+
  ylim(-1.5,6.5)+
  xlab("exp. de vida")+
  scale_color_manual(values = colores_desarrollo)

z2005<-ggplot(datos_knn_2005, aes(exp_de_vida, gdp, color=nivel))+
  geom_point()+
  ggtitle("Exp. de vida vs GDP, año 2005")+
  xlim(-3.2,2.2)+
  ylim(-1.5,6.5)+
  xlab("exp. de vida")+
  scale_color_manual(values = colores_desarrollo)

z2010<-ggplot(datos_knn_2010, aes(exp_de_vida, gdp, color=nivel))+
  geom_point()+
  ggtitle("Exp. de vida vs GDP, año 2010")+
  xlim(-3.2,2.2)+
  ylim(-1.5,6.5)+
  xlab("exp. de vida")+
  scale_color_manual(values = colores_desarrollo)

z2015<-ggplot(datos_knn_2015, aes(exp_de_vida, gdp, color=nivel))+
  geom_point()+
  ggtitle("Exp. de vida vs GDP, año 2015")+
  xlim(-3.2,2.2)+
  ylim(-1.5,6.5)+
  xlab("exp. de vida")+
  scale_color_manual(values = colores_desarrollo)

grid.arrange(z2000,z2005,z2010,z2015, ncol=2, nrow=2)

train =  datos_knn_2000[, -c(1,2)] # saco la columna con la etiqueta
test =  datos_knn_2000[, -c(1,2)]

labels_train = datos_knn_2000[,2] # guardo las etiquetas
clasif_datos <- knn(train = train, test = test, cl = labels_train, k = 3)
summary(clasif_datos)

##  Developed Developing 
##         30        126

labels_test = datos_knn_2000[,2]
confusionMatrix(as.factor(clasif_datos), as.factor(labels_test))

## Confusion Matrix and Statistics
## 
##             Reference
## Prediction   Developed Developing
##   Developed         24          6
##   Developing         6        120
##                                           
##                Accuracy : 0.9231          
##                  95% CI : (0.8695, 0.9596)
##     No Information Rate : 0.8077          
##     P-Value [Acc > NIR] : 4.742e-05       
##                                           
##                   Kappa : 0.7524          
##                                           
##  Mcnemar's Test P-Value : 1               
##                                           
##             Sensitivity : 0.8000          
##             Specificity : 0.9524          
##          Pos Pred Value : 0.8000          
##          Neg Pred Value : 0.9524          
##              Prevalence : 0.1923          
##          Detection Rate : 0.1538          
##    Detection Prevalence : 0.1923          
##       Balanced Accuracy : 0.8762          
##                                           
##        'Positive' Class : Developed       
## 

# knn original -----------------------------------------------------------------

#set.seed(73)

datos_knn_1 <- datos_knn


table(datos_knn_1$nivel)

## 
##  Developed Developing 
##        480       2049

barplot(table(datos_knn_1$nivel)/2559, col=c("deepskyblue3","firebrick"))

Nrep = 20
Ntest = 30
max_val = 20
valores = 1:max_val

resultados_crossval_train = matrix(NA, length(valores), 51)
resultados_crossval_test = matrix(NA, length(valores), 51)

datos_knn_1 <- datos_knn_1[,-1]


for(k in 1:max_val){
  
  for(i in 1:50){
    
    indices = sample(seq(1, nrow(datos_knn_1)), 0.2*nrow(datos_knn_1))
    test = datos_knn_1[indices, -1]
    train = datos_knn_1[-indices, -1]
    test_labels = datos_knn_1[indices, 1]
    train_labels = datos_knn_1[-indices, 1]
    
    clasificador_datos_k_train <- knn(train = train, test = train, cl = train_labels, k = k)
    clasificador_datos_k_test <- knn(train = train, test = test, cl = train_labels, k = k)
    
    exactitud_train = confusionMatrix(as.factor(clasificador_datos_k_train), as.factor(train_labels))$overall[[1]]
    exactitud_test = confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$overall[[1]]
    
    resultados_crossval_train[k, i] = exactitud_train
    resultados_crossval_test[k, i] = exactitud_test
  }
  
  resultados_crossval_train[k,51] <- mean(resultados_crossval_train[k,1:50])
  resultados_crossval_test[k,51] <- mean(resultados_crossval_test[k,1:50])
  
  
  if(k > 1){
    if(resultados_crossval_test[k, 51] > mayor){
      mayor <- resultados_crossval_test[k, 51] 
      matriz_confusion = data.frame(confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$table)
      titulo <- k
    }
  }else{
    matriz_confusion = data.frame(confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$table)
    titulo <- k
    mayor <- resultados_crossval_test[k, 51] 
  }
  
  
}


exact_train = resultados_crossval_train[,51]
exact_test = resultados_crossval_test[,51]

plot(valores, exact_train, col = 'darkblue', pch = 19, ylim = c(0.85,1), cex = 0.5, xlab = 'Valor de k', ylab = 'Exactitud', main = 'Exactitud en función del k')
lines(valores, exact_train, col = 'darkblue')
points(valores, exact_test, col = 'darkorange', pch = 19, cex = 0.5)
lines(valores, exact_test, col = 'darkorange')
legend('topright', legend = c('train', 'test'), col = c('darkblue', 'darkorange'), pch = 19)
abline(h = 0.9, col = "red")

plotTable <- matriz_confusion %>%
  mutate(goodbad = ifelse(matriz_confusion$Prediction == matriz_confusion$Reference, "good", "bad")) %>%
  group_by(Reference) %>%
  mutate(prop = Freq/sum(Freq))

#    fill alpha relative to sensitivity/specificity by proportional
# outcomes within reference groups (see dplyr code above as well as
# original confusion matrix for comparison)

ggplot(data = plotTable, mapping = aes(x = Reference, y = Prediction, fill = goodbad, alpha = prop)) +
  geom_tile() +
  ggtitle(as.character(titulo)) +
  geom_text(aes(label = Freq), vjust = .5, fontface  = "bold", alpha = 1) +
  scale_fill_manual(values = c(good = "green", bad = "red")) +
  theme_bw() +
  xlim(rev(levels(matriz_confusion$Reference)))

# Knn con menos datos de developing --------------------------------------------

datos_knn_2 <- datos_knn

table(datos_knn_2$nivel)

## 
##  Developed Developing 
##        480       2049

developed <- datos_knn_2[is.element(datos_knn_2$nivel,"Developed"),]
developing_0 <- datos_knn_2[is.element(datos_knn_2$nivel,"Developing"),]
developing <- sample_n(developing_0, 1000)
datos_knn_2 = rbind(developed,developing)
barplot(table(datos_knn_2$nivel)/1480, col=c("deepskyblue3","firebrick"))

Nrep = 20
Ntest = 30
max_val = 20
valores = 1:max_val


resultados_crossval_train = matrix(NA, length(valores), 51)
resultados_crossval_test = matrix(NA, length(valores), 51)

datos_knn_2 <- datos_knn_2[,-1]


for(k in 1:max_val){
  
  for(i in 1:50){
    
    indices = sample(seq(1, nrow(datos_knn_2)), 0.2*nrow(datos_knn_2))
    test = datos_knn_2[indices, -1]
    train = datos_knn_2[-indices, -1]
    test_labels = datos_knn_2[indices, 1]
    train_labels = datos_knn_2[-indices, 1]
    
    clasificador_datos_k_train <- knn(train = train, test = train, cl = train_labels, k = k)
    clasificador_datos_k_test <- knn(train = train, test = test, cl = train_labels, k = k)
    
    exactitud_train = confusionMatrix(as.factor(clasificador_datos_k_train), as.factor(train_labels))$overall[[1]]
    exactitud_test = confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$overall[[1]]
    
    resultados_crossval_train[k, i] = exactitud_train
    resultados_crossval_test[k, i] = exactitud_test
  }
  
  resultados_crossval_train[k,51] <- mean(resultados_crossval_train[k,1:50])
  resultados_crossval_test[k,51] <- mean(resultados_crossval_test[k,1:50])
  
  
  if(k > 1){
    if(resultados_crossval_test[k, 51] > mayor){
      mayor <- resultados_crossval_test[k, 51] 
      matriz_confusion = data.frame(confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$table)
      titulo <- k
    }
  }else{
    matriz_confusion = data.frame(confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$table)
    titulo <- k
    mayor <- resultados_crossval_test[k, 51] 
  }
  
  
}


exact_train = resultados_crossval_train[,51]
exact_test = resultados_crossval_test[,51]

plot(valores, exact_train, col = 'darkblue', pch = 19, ylim = c(0.85,1), cex = 0.5, xlab = 'Valor de k', ylab = 'Exactitud', main = 'Exactitud en función del k')
lines(valores, exact_train, col = 'darkblue')
points(valores, exact_test, col = 'darkorange', pch = 19, cex = 0.5)
lines(valores, exact_test, col = 'darkorange')
legend('topright', legend = c('train', 'test'), col = c('darkblue', 'darkorange'), pch = 19)
abline(h = 0.9, col = "red")

plotTable <- matriz_confusion %>%
  mutate(goodbad = ifelse(matriz_confusion$Prediction == matriz_confusion$Reference, "good", "bad")) %>%
  group_by(Reference) %>%
  mutate(prop = Freq/sum(Freq))

#    fill alpha relative to sensitivity/specificity by proportional
# outcomes within reference groups (see dplyr code above as well as
# original confusion matrix for comparison)

ggplot(data = plotTable, mapping = aes(x = Reference, y = Prediction, fill = goodbad, alpha = prop)) +
  geom_tile() +
  ggtitle(as.character(titulo)) +
  geom_text(aes(label = Freq), vjust = .5, fontface  = "bold", alpha = 1) +
  scale_fill_manual(values = c(good = "green", bad = "red")) +
  theme_bw() +
  xlim(rev(levels(matriz_confusion$Reference)))

# Knn con MUCHOS menos datos de developing -------------------------------------

datos_knn_3 <- datos_knn

table(datos_knn_3$nivel)

## 
##  Developed Developing 
##        480       2049

developed <- datos_knn_3[is.element(datos_knn_3$nivel,"Developed"),]
developing_0 <- datos_knn_3[is.element(datos_knn_3$nivel,"Developing"),]
developing <- sample_n(developing_0, 480)
datos_knn_3 = rbind(developed,developing)
barplot(table((datos_knn_3$nivel))/960, col=c("deepskyblue3","firebrick"))

Nrep = 20
Ntest = 30
max_val = 20
valores = 1:max_val

resultados_crossval_train = matrix(NA, length(valores), 51)
resultados_crossval_test = matrix(NA, length(valores), 51)

datos_knn_3 <- datos_knn_3[,-1]


for(k in 1:max_val){
  
  for(i in 1:50){
    
    indices = sample(seq(1, nrow(datos_knn_3)), 0.2*nrow(datos_knn_3))
    test = datos_knn_3[indices, -1]
    train = datos_knn_3[-indices, -1]
    test_labels = datos_knn_3[indices, 1]
    train_labels = datos_knn_3[-indices, 1]
    
    clasificador_datos_k_train <- knn(train = train, test = train, cl = train_labels, k = k)
    clasificador_datos_k_test <- knn(train = train, test = test, cl = train_labels, k = k)
    
    exactitud_train = confusionMatrix(as.factor(clasificador_datos_k_train), as.factor(train_labels))$overall[[1]]
    exactitud_test = confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$overall[[1]]
    
    resultados_crossval_train[k, i] = exactitud_train
    resultados_crossval_test[k, i] = exactitud_test
  }
  
  resultados_crossval_train[k,51] <- mean(resultados_crossval_train[k,1:50])
  resultados_crossval_test[k,51] <- mean(resultados_crossval_test[k,1:50])
  
  
  if(k > 1){
    if(resultados_crossval_test[k, 51] > mayor){
      mayor <- resultados_crossval_test[k, 51] 
      matriz_confusion = data.frame(confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$table)
      titulo <- k
      test$predichos <- clasificador_datos_k_test
      graf <- ggplot(test, aes(exp_de_vida, gdp, color = predichos))+
        geom_point()+
        xlab("exp. de vida")+
        scale_color_manual(values = colores_desarrollo)
      
      test2<-test
      test2$nivel <- test_labels
      graf2 <- ggplot(test2, aes(exp_de_vida, gdp, color = nivel))+
        geom_point()+
        xlab("exp. de vida")+
        scale_color_manual(values = colores_desarrollo)
  
        
        
    }
  }else{
    matriz_confusion = data.frame(confusionMatrix(as.factor(clasificador_datos_k_test), as.factor(test_labels))$table)
    titulo <- k
    mayor <- resultados_crossval_test[k, 51] 
    test$predichos <- clasificador_datos_k_test
    graf <- ggplot(test, aes(exp_de_vida, gdp, color = predichos))+
      geom_point()+
      xlab("exp. de vida")+
      scale_color_manual(values = colores_desarrollo)
    
    test2<-test
    test2$nivel <- test_labels
    graf2 <- ggplot(test2, aes(exp_de_vida, gdp, color = nivel))+
      geom_point()+
      xlab("exp. de vida")+
      scale_color_manual(values = colores_desarrollo)
  }
  
  
}


exact_train = resultados_crossval_train[,51]
exact_test = resultados_crossval_test[,51]

plot(valores, exact_train, col = 'darkblue', pch = 19, ylim = c(0.85,1), cex = 0.5, xlab = 'Valor de k', ylab = 'Exactitud', main = 'Exactitud en función del k')
lines(valores, exact_train, col = 'darkblue')
points(valores, exact_test, col = 'darkorange', pch = 19, cex = 0.5)
lines(valores, exact_test, col = 'darkorange')
legend('topright', legend = c('train', 'test'), col = c('darkblue', 'darkorange'), pch = 19)
abline(h = 0.9, col = "red")

plotTable <- matriz_confusion %>%
  mutate(goodbad = ifelse(matriz_confusion$Prediction == matriz_confusion$Reference, "good", "bad")) %>%
  group_by(Reference) %>%
  mutate(prop = Freq/sum(Freq))

#    fill alpha relative to sensitivity/specificity by proportional
# outcomes within reference groups (see dplyr code above as well as
# original confusion matrix for comparison)

ggplot(data = plotTable, mapping = aes(x = Reference, y = Prediction, fill = goodbad, alpha = prop)) +
  geom_tile() +
  ggtitle(as.character(titulo)) +
  geom_text(aes(label = Freq), vjust = .5, fontface  = "bold", alpha = 1) +
  scale_fill_manual(values = c(good = "green", bad = "red")) +
  theme_bw() +
  xlim(rev(levels(matriz_confusion$Reference)))

grid.arrange(graf,graf2,nrow=1,ncol=2)